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ABSTRACT 

The Burgers equation with a small viscosity term, initial and periodic boundary condi- 
tions is resolved numerically using a spatial approximation constructed from an orthonormal 
basis of wavelets. 

The algorithm is directly derived from the notions of multiresolution analysis and tree 
algorithms. Before the numerical algorithm is described these notions are first recalled. 
The method uses extensively the localization properties of the wavelets in the physical and 
Fourier spaces. Moreover, we take advantage of the fact that the involved linear operators 
have constant coefficients. Finally, the algorithm can be considered as a time marching 
version of the tree algorithm. 

The most important point is that an adaptative version of the algorithm exists: it allows 
one to reduce in a significant way the number of degrees of freedom required for a good 
computation of the solution. 

Numerical results and description of the different elements of the algorithm are provided 
in combination with different mathematical comments on the method and some comparison 
with more classical numerical algorithms. 


1 Research was supported in part by the National Aeronautics and Space Administration under NASA 
Contract No. NAS1-18605 while the author was in residence at the Institute for Computer Applications 
in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665 and also by 
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1. INTRODUCTION 


The nonlinear parabolic equation 

du udu vd 2 u 

di ^ dx dx 2 

known as the Burgers equation is one of the simplest combining both nonlinear propagation 
and diffusive effects. It represents a first step in the hierarchy of approximations of the 
Navier-Stokes equation. Solutions of this equation exhibit a delicate balance between the 
nonlinear advection and the diffusion terms. This situation leads to solutions with rapid 
and localized variations. Moreover, exact solutions are known (thanks to the Hopf-Cole 
transformation). Thus, this equation appears to be a very convenient and useful test problem 
for new numerical schemes. 

This paper is devoted to the numerical resolution of the Burgers equation (1) with pe- 
riodic boundary conditions on [0,1] and known initial condition: tz(0,i) = u(l, t); u(x, 0) = 
sin(27rx). We use a new algorithm based on a spatial approximation provided by an orthonor- 
mal wavelet basis exploiting as much as possible the numerical localization and regularity 
properties of the wavelets. One of the greatest advantages of this algorithm is that the 
whole method can automatically adapt itself to the generation of large gradient regions in 
the solution. 

After a short review of the notions of multiresolution analysis and tree algorithms that 
are required for the construction of orthonormal wavelet bases and that are the keystones 
of our numerical algorithm, some properties of r — regular wavelet bases and particularly 
of spline wavelets (Section 2) are recalled. The choices and definitions for the numerical 
algorithm are discussed first in the regular case (Section 3.1) and secondly in the adapted 
case (Section 3.2). Numerical results are presented in Section 4 with a complete description 
of the different elements of the method. Comparison with a classical Fourier pseudospectral 
method is provided. Section 5 is devoted to concluding remarks. 

2. CONSTRUCTION AND BASIC PROPERTIES OF ORTHONORMAL 
WAVELET BASIS 

The natural framework in which to construct bases of wavelets is given by the multireso- 
lution analysis. This new concept, elaborated by S. Mallat and Y. Meyer, is also well adapted 
to the approximation of functions for numerical purposes. 

The existence of a “Fast Wavelet Transform,” based on the class of tree algorithms studied 
by S. Mallat and I. Daubeschies, is also a nodal point. 
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We briefly describe the main definitions and results we will need. The reader is referred 
to the book of Y. Meyer [Meyer 1987] for a complete account on the subject. Local analysis 
properties will demonstrate the quality of the wavelet approximation and help to explain 
the flagging process used in the adaptative version of the algorithm. The periodic spline 
wavelets will be presented at the end of this section and the different properties motivating 
their choice for the numerical resolution of partial differential equations will be recalled. 

2.1. Multiresolution Analysis 


As introduced by S. Mallat and Y. Meyer, a multiresolution of L 2 (M n ) is a specific 
approximation scheme for finite energy functions. In this paper we restrict ourselves to 
n = 1. A multiresolution of L 2 (1R) is an increasing sequence of closed linear subspaces 
Vj,j E 2/ such that: 

+ oo -foo 

(2.1) Vj = {0}, (J Vj is dense in L 2 (1R ) 

— OO — OO 

(2.2) /(*) g Vj «=> /( 2x) e v j+1 

(2.3) f(x ) G Vo <=> f(x - Jfc) G Vo, VkEZ 

(2.4) There exists a function g(x) in Vo such that the collection g(x — k), k G ZZ is a Riez 
basis for Vo 

(A collection of vectors ej,j E J , in a Hilbert space H is a Riez basis if any vector x E H 
can be written in a unique way as a sum x = Y, a j e j where (53 \a.j\ 2 ) l l 2 is finite and defines 
an equi valient norm.) 

The multiresolution analysis is called r-regular if the function g(x) defined by (2.4) has 
the following property: 

(2.5) |d“s(x)| < C p (l + [x |) — p for a < r, all x E R and all integer p. 

This property can be interpreted as a quantification of the localization properties of the 
function g in the physical space (a = 0) and in the Fourier space (a > 1). 

Let us give two examples of multiresolution analysis that generate classical spaces. The 
first one, to which people familiar with spectral methods automatically think, is provided by 
the Paley- Wiener analysis. In this analysis, Vj is made of functions whose spectrum belongs 
to the interval [— 2 J_1 , 2 J_1 ]. g is the cardinal sinus g(x) = . It is well-known that g 
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is infinitely derivable but is badly localized in the physical space. Consequently, for every 
value of r, this multi-resolution analysis is not r-regular. 

The second example is the spline multiresolution which will provide the spline wavelets 
that will be used in the following numerical implementation of the algorithm . Here, given an 
integer m > 2, Vj consists of all functions with m— 2 continuous derivatives whose restrictions 
to any interval [&2 -J , ( k -f 1 )2 — J ] coincide with a polynomial of degree less or equal to m — 1. 
g{x) is given by the so-called basic spline g m (x) which is the m fold convolution of the 
characteristic function of [0, 1]. g m is supported by [0,m] and, satisfies (2.5) for r = m — 2. 

2.2. The Construction of Wavelets from the Multiresolution Analysis 

The construction of the wavelet basis stems from the fact that during the process of 
refinement in the approximation one wants to only store the improvement from the approx- 
imation j to the approximation j + 1. Mathematically, one introduces at each step j, the 
subspace Wj, defined as the orthogonal complement of Vj in Vj +j. The Wj space family 
satisfies the scaling (2.2) and translation invariance (2.3) properties imposed on Vj, so that 
attention can be focused on Wq. Then one has the fundamental theorem proved by S. Mallat 
and Y. Meyer. 

Theorem: There exists a function ip of Wq such that {ip(x — k), k £ 2Z\ is an orthonormal 
basis of Wo- ip has the same regularity properties (2.5) as g, and ip is an oscillating function: 

(2.6) J x k ip(x)dx = 0 for 0 < k < r. 

Then, the family of functions {^^(x) = 2^ 2 ip{Vx — k),k £ Z} is an orthonormal family 
of Wj and, as L 2 (IR ) = OjezrWj, {ipjk(x),j £ 2Z,k £ is a hilbertian basis of L 2 (1R). 
This family is an orthonormal wavelet family. Combined with (2.5) for a > 0, equation (2.6) 
characterizes the localisation of the wavelet generating function in the Fourier space. 


The proof of the theorem introduces an orthonormal basis of Vj that will be denoted 
by {<pjk(x) = 2?l 2 <p{2?x — k),k £ %}. These functions are called by S. Mallat the scaling 
functions due to the fact that / <p(x)dx = 1. They satisfy the same regularity properties (2.5) 


as g. 


Given an integer j 0 and writing L 2 (JR ) = Vj 0 © Y,j>j 0 W) one obtains another hilbertian 
basis of L 2 (M ) : {<pj 0 ,k,i>ik,j > jo,k £ Z}. 

Let us go back now to the two examples: 

The Paley-Wiener Analysis provides the wavelet generating function ip(x) given by: 


ip(x) = 


2 cos 7rx — sin 2irx 
7 r 2x — 1 


and ip(u ;) = e ‘™Xi/ 2 <|u,l<i 
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where ip stands for the Fourier transform of ip : ip(w) = / ip{x)e~ 2xirwx dx. 

ip has a compact support but the localization of ip is poor as it decreases as ^ at infinity 
and ip(uj) — 0 if |u/| < which means, in a weak sense, that all the moments of ip are 
vanishing. 

The spline multiresolution analysis provides the “exponentially localized wavelets” de- 
rived by P. G. Lemari4 and G. Battle [Battle 1987, Lemarie 1989]. The generating wavelet 
ip(x) satisfies the following properties: 


(2.7) 


ip{x) and all its derivatives to the order m — 2 axe continuous 


(2.8) 


there exist e > 0 such that \ip(x)\ < Coe e ' x ', l^~( x )l — ^i e • • ' » 

l^r 1 ' * < 7 -’ e "* W 


(2.9) 


/ +oo 

x k ip(x)dx = 0 for k — 0, ■ ■ • , m — 1 

-OO 


ip is given in the Fourier space by 


Tp{w) 


rim - 1 (°°g ?) 


(=)" l p — i ( sin ’ (=?)) 

where P m is the m order polynomial given by 


sm 


2m 


/7TW\ 

vT7 


Pm- i(sin 2 u;) _ y, 1 
(sin w) 2m (w + fcir) 2TO 

The scaling function <p is also “exponentially localized” and satisfies: 

« 1 sin m (7riw) 

W ~ (™) m [P m _i(sin 2 (7rw))]3' 

For m = 6 the generating spline wavelet ip and scaling function <p are plotted in Figure 1. 


2.3. Tree Algorithms 

The tree or pyramidal algorithms are the basic tools for the fast computation of wavelet 
coefficients. Also, the numerical algorithm presented here is efficient thanks to its pyramidal 
structure. 
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If jb f is a given integer, fixing the level of approximation at 2 -jM , the tree algorithm 
allows fast computation of the wavelet coefficients from jm — 1 to an order jo , jo < j^, 
starting from the scaling coefficients in Vj M . 

Writing 

JM~ 1 

(2.10) V ju = V K © Y. W i 

j-jo 

and starting from the j\i approximation of a function U, we have 


Hvj„(C0(x) = E - k). 

k£,Zs 

(Generally, II x stands for the orthogonal projector from L 2 (M ) on the subspace X.) At 
each level j, S. Mallat [Mallat 1988] shows that, due to the orthogonality properties and the 
scaling, TLyp.^iJJ) and IIv' J _ l ({7), written as: 

= E ‘ ) 

lez 

and 

nr,.,(tf)(s) = E 

can be computed from II Vj{U) using the formulas: 


dj-itl — 9$£ k ) 

k 


and 

= J2 Cj,k h(2t - k) 
k 

where g and h are two discrete functions depending only on the multiresolution analysis and 
independent of jf. These functions are precisely defined as: 

Vn <E N : g{n) = -^= < V’-i.o, 4>o,n > and h{n) = -^= < <f> 0 , n > 


with 


5(~ n ) = a{ n ) ) h(-n) = h(n) 


where < ., . > stands for the L 2 (1R) inner product (< /, g >= / f(x)g(x)dx). 

Practically, it means that with the storage of two discrete functions h and g which are well 
localized thanks to the localization of ip and <p , one computes with an order of NLogN(2 jM = 
N ) operations the (2 JM — 2 70 ) wavelet coefficients and the 2 J0 scaling coefficients corresponding 
to the formula (2.10). The discrete functions h(n ) and g(n) are plotted on Figure 2 in the 
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physical and Fourier spaces in the case of the spline multiresolution with m = 6. They 
classically appear as a low pass filter and a high pass filter respectively. 

2.4. Approximation properties of Wavelets 

Approximation properties of orthonormal wavelets can be interpreted from different 
points of view: 

As the space generated by the wavelets up to an order j — 1 is the space Vj of the 
multiresolution, the approximation results are driven by the properties of the scaling function 
2^ 2 <f>(2^x) that, with its translates 4>jk, generate Vj. It is then a more classical problem of 
approximation by translates. The general results relating the quality of the approximation 
to the degree r of r-regularity of the multiresolution can be found in Y. Meyer’s book [Meyer 
1990]. 

Moreover, if the space Vj is split in the different contributions of W), / < j — 1, a char- 
acterization of the Sobolev space H‘(JR ) can be obtained from the decay of the wavelet 
coefficients as : 

Theorem: If f belongs to H~ r (IR), ifVj,j G Z is a r-multiresolution of L 2 (M) and if —r < 
s <r then f belongs to H t (M) if and only i/IIy 0 (/) G L 2 (M) and ||IIty,(/)||2 = £ IN, 

where £j G l 2 (W). 

Another important result for approximation by wavelets deals with the characterization 
of singularities from the wavelet decomposition. The singularity of a function at a point x 0 
can be described with the Lipschitz exponent a defined as: 

f(x) is a Lipschitz at x 0 ,(0 < a < 1) if and only if for all x in a neighborhood of x 0 one 
has: 

\f(x)-f(x 0 )\ = 0(||® -XolT). 

The following theorem holds: 

Theorem: A function f belonging to L 2 (M) is a Lipschitz with 0 < a < 1 in all points of 
an open interval if and only if for all x in the interval one has: 

I < Mil. > I = 0(2-'<“ +1 «). 


As pointed out by Stephane Mallat [Mallat 1990], an extension of this result can be 
obtained when the initially a singular function f is smoothed by a smoothing kernel. For 
instance with a gauusian of variance a one has the following result: 
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The wavelet coefficient of the smoothed function at the scale 2 * is of the same order of 
magnitude as the wavelet coefficient of the original function f at the scale s = v2 3 . 

The previous theorem then holds by changing 2 -J in \/2 -2j + cr 2 . 

On one hand, when 2 — ^ c, i.e., when the scales under study are larger than the 

smoothing scale, the previous theorem shows that the characterization of the degree a from 
the wavelet coefficients of the smoothed function is still possible. On the other hand, when 
2~ J <c a, the scaling effect disappears and no characterization of the degree of singularity is 
possible from the wavelet coefficients of the smoothed function. 

2.5. Periodic Spline Wavelets 

As pointed out by Y. Meyer [Meyer 1990], the complete tool box built in L 2 (1R) can be 
used in the periodic case for T 2 ([0, 1]) by introducing a standard periodization technique. It 
has been implemented by V. Perrier and C. Basdevant [Perrier and Basdevant 1989]. This 
technique consists at each scale in folding, around the integer values, the wavelets V'ifc and 
the scaling functions <f>j k centered in [0, 1] resulting in 

V’PjkC®) = J2 “ 0 

lez 

and 

4>vA x ) = ^;*( x - 0- 

lez 

This folding provides a multiresolution analysis of L 2 ([0, 1]). Due to the finite size of [0,1], 
j takes only positive values and k takes a finite number of values at each scale. 0 < k < 
2 J — 1. Moreover, the periodized function 4>po,o that generates Vpo is constant and equal to 
1. Following V. Perrier and C. Basdevant, the periodic m order spline scaling functions and 
wavelets are „ „ , , 

' oM = M (?) 

i 

„ V’Pj.oM = wit (f) 

where <f>pjk and V'Pjit generate respectively V p j and W p j. 

For simplicity in the notation, from now on the subscript p will be removed. When 
necessary the subscript up will be added for the non periodized functions. 

As far as numerical application is concerned, the spline multiresolution offers the following 
advantages. 

First, it is known that if the degree m of the spline is even, there exists a Lagrangian 
interpolant function gj M of Vj M such that gjM{k2~* M ) = £fc,o- (£ stands for the Kronecker 
symbol: 8 k>0 = 1 if and only if jfe = 0.) The functions g ]M {x - 0 <k< are then 
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an unconditional basis of V*- Thanks to this basis, it is possible, for any function U(x), to 
consider its interpolation on Vj M given by 

n.»i M (tO(.)= E v(*2-*'),m,(x - «-<“). 

0<fc<2 , Af 

This collocation projection is very useful when one wants to apply the classical but efficient 
pseudospectral treatment to nonlinear terms (see Section 3.1.5). 

A second point of interest stems the numerical localization of the generating functions 
V’up and <f> up in physical and Fourier spaces. Numerically, the exponential decay of these 
functions can be described by the required support to which V’up or <f> up (or ip up and <£ up ) can 
be restricted without altering their L 2 norms (equal to 1) up to a fixed precision a. 

In the physical space, the support of the periodized functions <f>jk and V’* is directly 
connected to the support of <f> up jk and V’up jk as soon as the length of these last supports is 
less than 1 (otherwise the support of <f> jk and ip jk is [0,1]). Moreover if S^ p and Sy, up are 
the supports of <f> up and V’up the supports of <f> up jk and V’up * are respectively and ^a. 
For m = 6, the following estimations are obtained: 


precision a 

support of <f> up 

support of V'up 

10 -3 


x] < 3.50 

|x| < 3.25 

1 

0 

H 

x\ < 4.75 

x < 4.75 


and an estimate of the increase of the length of these supports with m can be found in [Perrier 
and Basdevant 1989]. For j > 3,m = 6 and up to a precision a = 10 -4 , the functions (<f> up )jk 
and (V’up)* have a numerically compact support of length less than 1. 

In the Fourier space, the periodization does not affect the support of the functions p and 

A A A 

V’- If Si are the bounds of the support of <f> (or ip), 2 j Si are the bounds of the support of <j> ik 
(or V»*). For m = 6 we obtained 


precision a 

support of <j> 

support of V> 

HP 

w <0.625 

0.34 < |u/| < 1.29 

h-* 

O 

1 

|to| < 0.685 

0.27 < M < 1.40 


The length of these supports decreases as m increases. We will strongly use the fact that, 
for m = 6 and a = 10" 3 , the numerical supports of ^j k and i’j+ 2 ,k are disconnected. 

A third advantage of the spline wavelets is the very easy control of their cancellation 
(the first m — 2 moments of each wavelet vanish). As has been shown in Section 2.4, on 
one hand, this property is required to have a good convergence of the wavelet series and a 
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good characterization of the local properties of the analyzed function. On the other hand, 
as m increases the localization in the physical space of the wavelet decreases. To obtain a 
“satisfactorily localized function” with “good cancellation,” one has to define a compromise 
for the choice of m . All the numerical results presented in this paper are obtained with 
m = 6 . 

To be complete, one must recall that an efficient implementation of periodic spline wavelet 
decomposition has been described by V. Perrier and C. Basdevant [Perrier et al. 1989]. It 
allows computation with order of N Logi^N') operations of the wavelet coefficients of a periodic 
function described by N regularly distributed collocation points. This algorithm has been 
extended to a non complete (non regular) distribution of points. 

3. THE ALGORITHM 

Taking apart the Lagrangian methods, numerical algorithms can be roughly speaking 
separated into three families: finite difference schemes, finite element methods and spectral 
schemes. 

The class of finite difference schemes is known to be very efficient in intricate configura- 
tions because it can be easily adapted to resolve localized difficulties such as large gradient 
regions or sharp boundary conditions. This flexibility, however, hides a difficulty in control- 
ling precision. Moreover, the computation cost increases considerably with the precision. 

The finite element schemes offer a better controlled precision and are extensively used 
for intricate boundary condition problems in multidimension. However they suffer from 
difficulties in getting very flexible and fast algorithms. 

The precision of spectral methods is well-known ([Gottlieb and Orszag 1977], [Canuto 
et al. 1987]) and their efficiency for problems with a high degree of regularity has been 
demonstrated. However, they are fundamentally ill-suited to problems that develop strong 
local gradients or discontinuities in their solutions. 

The wavelet approach provides a multiscale decomposition based on orthonormal, regular 
and numerically well localized functions. It can then a priori offer an interesting compromise 
between precision, efficiency, and adaptability. 

3.1. Description of the Algorithm on a Regular Grid 

3.1.1. The Grid Points 

Given a family of wavelets, one can characterize the space X that they generate by the 
dyadic grid built from the set of all y jk = £ + ^ such that ^ jk belongs to X (see Figure 3). 
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If X is a space Vj M , then a regular grid corresponding to the centers of the scaling functions 
is built from the X*. is the level of approximation. 

3.1.2. The Algorithm 


The algorithm is first described generally for the evolution equation 

' *L + LU + G{U) = 0 


(3.1) 


U = Uo for t = 0 


( t >0,0 <x < 1 


where L is a linear operator and G is a function of U and its derivatives. An approximate 
solution U M , belonging to a trial space X M , (X M = Vj M in the case of a regular grid) is 
sought as a solution of the equation: 

j + L m Um = U Xm (G(U m )) 

1 U M (x } 0) = U jM U 0 . 

Here Lm is an approximation of L, of the form Px M LHx M , where IIy M is the orthogonal 
projection onto Xm, and Px M is another projection onto Xm parallel to the orthogonal of 
Ym, another space to IIx Ar (G r (uAf)); the approximate equation becomes 


(3.2) 


' ^ + LuUm = G m 

, 0) = t^oj^a;). 


Using the classical name conventions for the Methods of Weighted Residuals (MWR) (see 
[Canuto et al. 1987]), we call trial functions the wavelets V’jfcjO < j < j‘m,0 < k < 2 J and 
<f>o,o and introduce a family of functions 6 T] t G T, called the test functions: they generate our 
space Ym- The weighted residuals minimization statement is then equivalent to the following 
variational formulation: 


(3.3) 


< + LmUm — Gm, v >= 0 for all v G Ym 

k Um(x,0) = U Om (x) 


where < •, • > stands for the L 2 inner product on [0,1]. In the case of the Burgers equation 
(1), we will have LU = — v ^p' , G(£7) = — • * , Gm{U) = II C (G(U) and the test functions 
8 r will be defined in the next paragraph. 
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.1.3. Choice of the Test Functions and of the Time Discretization Scheme 

The general idea of the algorithm is to comply with the localization of the wavelets in 
the Fourier space and in the physical space. We know that, in comparison to the Fourier 
functions e tkx , the localization in the physical space implies the mixing of the frequencies 
that prevents the wavelets from being eigenfunctions of differential operators. However, they 
are not far from being so. Indeed, the localization of the wavelets in the Fourier space leaves 
the differential operators nearly diagonal. 

In the particular case of the spline wavelets, it has been shown previously that, to a given 
precision a, the numerical supports of ipjk and Vyi are disconnected as soon as | j —j'\ > n(a) 
(for a = 10~ 3 ,n = 1). If D is a differential operator with constant coefficients, the same 
property holds also for Dipjk and ifij'k with a new precision a 1 depending on the spectrum of 
D. If, for instance D = (1 — i/Jp-), the precision is unchanged. 

The algorithm we are going to describe here takes large advantage of this property. 
Moreover, thanks to the fact that the coefficients of the Burgers equation are constant, some 
precalculations can be performed once and for all. Indeed, the main part of the work to be 
performed at each time step of the MWR can be done only one time and reused at every 
time step. This is done by a proper choice of the test functions 0 T . 

Let us be more precise in the case of an Euler time discretization scheme, implicit for the 
dissipation term and explicit for the convection term. This choice is made only for sake of 
simplicity in the following presentation and more sophisticated schemes are currently used 
for the numerical implementations (see Section 4). 

Equation 3.3 is discretized as follows: 

' < (1 - >=< US, - > for all v € Yu 

i 

. U M {x, 0) = Uq m {x). 

Choosing r = (j. As) and 6 r = djk such that 
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that can be written as: 


(3.6) 


' < off',** >=< u&,e, k > + < (vz,)\o' ik > 

> < ~ oo >— < Ufa,6 0 > 


where we introduced another set of functions: 


(3.7) 


0, jk( x ) = -f (i ( 7 - vAt &)) 1 

. 0’ 0 {x) = 0. 


Using a Crank Nicholson discretization scheme and spline wavelets of order m = 6, the 
function 6 associated with ip up using (3.5) is plotted on Figure 5. Roughly speaking, the 
functions look like the corresponding wavelets ipjk'- In the case of the spline wavelets, the 
6jk functions are exponentially localized in both Fourier and physical spaces. Moreover they 
have the same number of zero moments as the original wavelets and their numerical support 
in the Fourier space is included in those of the corresponding wavelets. At each scale, the 
Ojk deduce one from the other by translation but the wavelet rescaling property from one 
scale to the other is no longer valid. 


3.1.4. Computation of < Ufa, 6j k > 

A fast computation of these coefficients is possible thanks to the previous remarks. If, 
following 2.10, one writes: 

vj < j'm - 2 u VjM = h Vh7 © X) n Wl 

j M -l>l>j+2 

following 3.1.1 with D = (1 — j/Af J^-), a numerically good approximation of < Ufa, 6jk > is 
<n v^(Ufa),e jk >. Then, recalling that: 

ny. +a = u Vj+l © n Wj+l 

one gets 

< Ufa, 6jk >^< n v i+2 (Ufa),e jk >~< Uv s+l (ufa),9 jk > + < n w . +l (ufa),e Jk > . 

Since we have Hv j+1 {Ufa) = E Cj+i.^i+U and Hw Hl {Ufa) = E d ]+hl tpj + u we obtain 

< Ufj i @jk >— k ~ 2^) + — 2k) 

L 

where ctj and fij are two families of discrete functions parameterized by j and depending on 
the multiresolution analysis, on the equation and on the time scheme. They are defined as 

Vn € N ocj(n) =< dj_ li0 > and j3j(n) =< d 3 _ lfi > 
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and 

a(—n ) = ay(n) and (3j(—n ) = (3j(n). 

As in the tree algorithm (Section 2.3), the calculation of {< Ufa, Qjk >, for 0 < j < 
3m-\ an d 0 <k < V} involves two filtering procedures. The whole calculation requires 
order JVXogjV operations. The coefficients of the filters ay and /3j are stored in the memory 
at the beginning of the calculation in an array of size proportionnal to NLogN . 

One wants to emphasize that the one-time computation involved in the calculation of the 
functions 6j k is possible because the operator I — vAt-^s has constant coefficients. 

The computation of < Ufa, 6'j k > can be done the same way using two other families of 
filters, ay-,/3', as soon as (Ufa) 2 is computed (see next section). 

Recalling that the filters ay , /?j, a!-, /?'• depend on At, it must be said that the independence 
of At with respect to j is not required. Practically it means that the time step At can be 
adapted to the spatial scale, a possibility that can be useful in optimizing the time stability 
and precision conditions. 

3.1.5. Estimation of the Nonlinear Terms 


The calculation of the nonlinear terms, precisely (Ufa) 2 , is done using the classical collo- 
cation technique. 

3.1.6. Implementation of the Algorithm 

At each time step, the algorithm works as follows: supposing that (Ufa) is known by its 
values at grid points, its scaling coefficients or its wavelets coefficients, Ufaf 1 is computed 
using an order of NLogN operations as shown in Table 1. 


u n „e 


v iu 


Pseudospectral 

► (U n „f 6 V iM 


♦ 


Tree Algorithm 


t 



> 
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3.2. Description of the Algorithm on a Non Regular Grid 


Up to now, Xm was equal to Vj M which means that every wavelet i/’jfc with(j, k) such 
that 0 < j < and 0 < k < V — 1 was considered in the approximation of the solution. 

However, the distribution of scales significantly building a given function in the sense 
provided by its wavelet decomposition has no reason to be regular in the whole interval 
[0,1]. Figure 7 shows the wavelet decomposition of the solution of (1) as large gradient 
regions develop locally around x = .5. According to the values of the wavelet coefficients, 
small scales ( j > 5) are only meaningful in a region very close to .5 and for time larger than 
.5/ir. More precisely, if we choose a given precision a on the L 2 norm and look at the number 
of wavelets required to approximate the solution at different times it appears that only a 
small number of coefficients of the wavelet decomposition are needed. For instance, with a 
L 2 norm precision of 10 -6 the number of wavelets required to reconstruct the solution at 
t = 1/x is 74. The distribution of these “active” wavelets follows the cone shape structure of 
Figure 7. Our main goal is to compute only the wavelet coefficients < UjfriV’jfc > in a region 
as close as possible to this “active” space. 

The adaptative version of our algorithm is designed in order to fulfill efficiently this 
requirement. It involves different steps: 

Step 1: Estimation of the required space of approximation for the next time step: In 
a way comparable to what is done in spectral methods (see Gottlieb - Orzag 1977), the 
accuracy of the computation is examined, locally in the case of wavelets, from the shape of 
the local spectrum provided by the wavelet decomposition at every point. 

Step 2: Adaptation of the space of approximation for the next time step: According 
to the result of the previous estimate and the knowledge of the operators of the Burgers 
equation, the space of approximation Xm is optimized (Xjf 4 ) for future time by truncation 
of non “active” small scale wavelets or local addition of smaller scale wavelets (that are 
supposed to become “active” during the future time steps). (The smallest scale wavelets in 
X Ma belong to Wj Ua - 1 ). 

Step 3: Time advancing on the adapted space Xm b - 

We are going now to describe in more detail the adaptative version of the algorithm 
in terms of flagging criterion (Step 1) and new trial and test functions (Steps 2 and 3). 
Moreover it will be shown that the general structure of the algorithm remains unchanged. 

3.2.1. New Trial Space X_m and Test Space Ym but Same Algorithm 

Figure 4 shows the selected values of ( j , k ) of an adapted grid generated by the adapted 
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version of the algorithm for the Burgers equation (1) at t — 1/7T. Precisely, 

X Ma = span {ip jk , (j, k) G JK} 
where 

JK = {(j, k ) such that 0 < j < <k <2* 
or such that j = jMi ki < k < k 2 
or, j = ;'jw + l,k[ <k< k' 3 } 

with Jm — 6, k\ = 22 = 31 k[ — 54 k 3 = 73. jm, = 8. 

To define completely the algorithm, one has also to define the ljw 0 space. This is done in 
exactly the same way as it has been presented in the regular case: the generating functions 
of Y Ma axe derived from the basis functions of Xmo using (3.5). Due to the local form of 
the algorithm, the whole procedure remains unchanged when the trial and test spaces are 
adapted. The complete procedure is identically reproduced without modification. 

This highly attractive property is due to the localization properties of the wavelets and 
of the algorithm in physical and Fourier spaces. 

The adaptive procedure is precisely a controlled modification of the space of approxima- 
tion in which one looks for the solution. As this procedure consists in adding or subtract- 
ing orthonormal functions, the so-called restriction and interpolation operators are exactly 
known: they are directly provided by the tree algorithm. Moreover, no reversibility problems 
are encountered during the restriction and interpolation procedures. 

For completeness, some words must be said on the flagging procedure used to modify 
the subscript space JK (this is the practical way to adapt Xm )• In the case of the Burgers 
equation, this procedure is very simple and cheap. Basically, starting from the smallest scale 
available (jss), the attention is focused on the smallest scale projection of the solution 

n w iss (Ub)= E < UM^jss.k > Vbss.*- 

k such that (jsStty^JK 

As shown in Section 2 and Figure 8, this function is localized in the strong gradient 
regions. The energy associated with the wavelets Vbssfcj precisely | < ipjss.k, > | 2 > 
reveals the amount of energy of the solution in the scale jss axound the point 2 iss+i 

The Burgers equation is such that, due to the presence of nonlinear terms in the 
equation that generate smaller and smaller scale, only small values of these coefficients are 
relevant (otherwise oscillations occur). Consequently, if | < i>j ss k, Ufa > | 2 is greater than 
a predetermined value, the space Xm of approximation is adapted by adding new smaller 
scaled wavelets ipj^i > C?a = jss + 1)> around the point 2 , S 5 +i + Adaptivity means also 
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restriction of the space x Um and, using a comparable criterion, small scale (jss) wavelets can 
be removed in a region where the terms | < Ufa, i>j ss t \ 2 are not significant. (This occurs in the 
second regime of the Burgers solution when the viscosity smooths the gradients previously 
generated during the first regime.) 

4. NUMERICAL RESULTS 

This section is devoted to a practical and numerical description of the algorithm. The 
basic elements of the problem and the algorithm are repeated for clarity. 

The periodic regularized Burgers equation 

r 8U . U8U _ v&U 

dt ' dx dx* 

t> 0,SG [0,1],!/ = 1211 

< 

U(0,t) = U(l,t) 

U(x,0) — sin(27rx) 

is solved numerically in the space of splines of order m - 6. The trial functions are the 
periodic spline wavelets i/> jk and 4> 0 o • The test functions are defined following the general 
procedure described in Section 3.1.5. All the following results are obtained using an Adams- 
Bashforth time scheme for the nonlinear convection term and a Crank- Nicholson time scheme 
for the diffusion term. The integration time step, for every scale is At = 10 -3 . We recall here, 
that the choice of the time scheme and the value of At directly influence the definition of the 
functions 8j k and However, knowing the time discretization scheme, their estimation is 
staightfor war d . 

The filters used in the wavelet decomposition are noted h and g. The filters defined to 
compute the scalar products < U^,6 jk > and < (U^) 2 ,6' jk > are noted a } ,(3 : and a',/3' 
respectively. The Figures 1, 2, and 5 show in the Fourier and in the physical spaces some 
selected elements of the different mathematical beings defined above. 

The evolution of the approximated solution U(x,t ) from t = 0 to t = ^ is plotted on 
Figure 6 where — 8 and a regular grid (Figure 3) is used. Considering the theoretical 
solution of the Burgers equation (see [Basdevant et al. 1986]), this can be considered as a 
reference estimation of the solution. Figure 7 shows the wavelet coefficients of the solution 
at the same times. Figure 8 shows the time evolution of the energy carried by the scale 

j - 6. 

At the same times, the same elements are plotted for a calculation run using jm — 

7 and a regular grid. Figure 9 reveals that the lack of resolution in the strong gradient 
region (around x = .5) induces numerical oscillations of the solutions close to those observed 
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when a Gibbs phenomenon occurs. These oscillations are however localized in space mainly 
thanks to the filtering operation performed by the spline projection of the variables and their 
derivatives. Figure 10 shows again the time evolution of an approximate solution obtained 
with a computation started with = 6 (that is to say 2 6 = 64 degrees of freedom and 
j ss = 5) and that has been obtained using the adaptative procedure. From a time close to 
t a — .5 fir when the local gradients begin to be very large the adaptative procedure adapts 
progressively the space to jM a = 7 and (ki,k 2 ) — (22,31) (that is to say 64 + 20 = 84 
degrees of freedom) and then (for .66/ 7r < t < 1.4/7t) to = 8 and (k[,k 2 ) = (54,73) 
(that is to say 64 + 20 + 20 = 104 functions). The corresponding dyadic grid is the one 
of Figure 4, the smallest scales j = 6 and j — 7 being used only when required. Finally 
on Figure 11, the evolution of the approximate solution using the same time discretization 
and a Fourier pseudospectral spatial discretization is plotted. The resolution is 2 7 = 128. 
First, one observes, in comparison to Figure 9 which is obtained with the same resolution 
(i.e the same number of degrees of freedom), that the oscillations are much more spread in 
the Fourier case, due to the non localization of the basis. Secondly, comparison with Figure 
10 shows that with a maximum of 104 well distributed degrees of freedom, a much better 
result can be obtain from the adapted algorithm. It is known (see [Basdevant et al. 1986]) 
that to improve the precision of a Fourier spectral method without increasing the number of 
degrees of freedom, intricate coordinate transformations are required that imply moreover 
an a priori knowledge of the solution shape. 

A comparison of the four different approximations presented above is made in Table 2. 
Tmn is the time needed to reach the maximum slope at x = 0.5, and S max is the value of 
this slope. No optimization study of the time marching scheme has been made that would 
have improved the estimates of iS'max an -d Tm** (see [Basdevant et al. 1986]). 

Table 2: Comparison of some different methods. 


Algorithm 

(AB-CN Time scheme) 

m 

SW 2. 

Exact: 152 .005 


Degrees of 
freedom 

Present Algorithm, m — 6 
3m = 8, regular grid 

1.64 

150.3 

No 

oscillations 

256 

Present Algorithm, m = 6 
3m = 7, regular grid 

1.63 

135.0 

Localized 

oscillations 

128 

Present Algorithm, m = 6 
3m = = 8 adapted grid 

1.64 

150.3 

No 

oscillations 

< 104 

Fourier pseudospectral 
n = 128 

1.62 

134.8 

Spread 

oscillations 

128 
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5. CONCLUSIONS 


The numerical resolution of the ID periodic and regularized Burgers equation has been 
performed using a new algorithm. This algorithm is based on the wavelet representation of 
the space of approximation. In the so-called regular grid case, this space is nothing more 
than the classical spline space, but in the non regular grid case it is more flexible than a 
classical grid refinement generated space. 

The special feature of this algorithm that uses explicitly the wavelet basis stands on 
the fact that it extensively takes into account the localization properties of the wavelets in 
the physical and Fourier spaces. Indeed, the exponential decay of the spline wavelets in 
both spaces allows, on one hand handling the differential operators as if they were nearly 
diagonal and, on the other hand, working in the physical space with very flexible orthogonal 
and numerically well localized functions. The scale decomposition provided by the wavelet 
coefficients and the knowledge of the Burgers operators are used to manage at each time 
step the space of approximation. The whole algorithm can then adapt itself to the solution. 
Moreover, in the case of constant regular operators, advantage can be taken by pre-inverting 
the problem at each scale once and for all at the beginning of the calculation. Then, the 
whole algorithm can be seen as a time advancing pyramidal algorithm. 

In this paper, the method has been described and a numerical application has been 
provided that uses spline wavelets of order 6. Some mathematical remarks have been given 
but the results on numerical analysis of the method will be published elsewhere. 

To date, we are not sure that the choice of the spline wavelets is optimal for efficient 
and adaptable resolution of partial differential equations; furthermore, the existence of a 
stronger connection between the wavelets and the resolved equation at each scale has not 
been investigated. However, it must be recalled that the spline wavelets achieve an optimal 
combined localization in physical and Fourier spaces, as fax as orthonormal wavelets are 
concerned. 

The limitation to a 1-D problem and periodic boundary conditions is not connected to 
the algorithm or the wavelet approach even if various specific problems will have to be faced. 
This choice has been made for clarity and simplicity. Extensions to more realistic problems 
involving different boundary conditions and multidimensions are currently underway. 
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Figure 2. a) h function, for the m = 6 spline multiresolution analysis, in 

physical and Fourier spaces. 

b) g function, for the m = 6 spline multiresolution analysis, in 
physical and Fourier spaces. 
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Figure 3. Regular dyadic grid (centers of rf> jk , 0 < j < j M ,0 < k < 2 j and of 

<t>o,o) {]M = 8 ). 
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Adapted dyadic grid j M = 6 ,{j = 6,*, = 22, k 2 = 31 }, {j = 7,k[ = 
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Figure 6. Time evolution of the approximated solution: m = 

regular grid algorithm. 256 degrees of freedom. 
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of the Approximated Solution 256 degrees of freedom 
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Figure 8. 


Time evolution of small scale (j = 6) energy. 
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Figure 10. 


Time evolution of the approximated solution m = 6 , ja / = 6 , j^a = 
7, adapted algorithm . The maximum number of degrees of 
freedom during the whole computation is 104. 
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